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ABSTRACT 

Observations of the high-mass star forming region AFGL 2591 reveal a large abundance of C0+, a 
molecule known to be enhanced by far UV (FUV) and X-ray irradiation. In chemical models assuming 
a spherically symmetric envelope, the volume of gas irradiated by protostellar FUV radiation is very 
small due to the high extinction by dust. The abundance of CO^ is thus underpredicted by orders of 
magnitude. In a more realistic model, FUV photons can escape through an outflow region and irradiate 
gas at the border to the envelope. Thus, we introduce the first 2D axi-symmetric chemical model of the 
envelope of a high-mass star forming region to explain the 00+ observations as a prototypical FUV 
tracer. The model assumes an axi-symmetric power-law density structure with a cavity due to the 
outflow. The local FUV flux is calculated by a Monte Carlo radiative transfer code taking scattering 
on dust into account. A grid of precalculated chemical abundances, introduced in the first part of this 
series of papers, is used to quickly interpolate chemical abundances. This approach allows to calculate 
the temperature structure of the FUV heated outflow walls self-consistently with the chemistry. 
Synthetic maps of the line flux are calculated using a raytracer code. Single-dish and interferometric 
observations are simulated and the model results are compared to published and new JCMT and 
SMA observations. The two-dimensional model of AFGL 2591 is able to reproduce the JCMT single- 
dish observations and also explains the non-detection by the SMA. We conclude that the observed 
CO^ line flux and its narrow width can be interpreted by emission from the warm and dense outflow 
walls irradiated by protostellar FUV radiation. 

Subject headings: stars: formation ~ stars: individual: AFGL 2591 - molecular processes - ISM: 
molecules - methods: data analysis 
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1. INTRODUCTION 

Forming high-mass O and B stars with a surface 
temperature larger than k, 10^ K emit most of their 
radiation in far UV (FUV) wavelengths. After a not 
yet well defined point of the evolution, also X-rays are 
emitted. When still deeply embedded in their natal 
cloud, most of this high-energy radiation is absorbed by 
the large column density of gas and dust towards the 
protostar. It is thus not available to direct observation 
but does influence the composition of the molecular 
envelope through heating and photoionization. An es- 
sential ingredient of star formation are bipolar outflows 
transporting the excess angular momentum outwards. 
When the fast, low-density gas of molecular outflows or 
jets expand into the surrouding molecular cloud, large 
cavities may result (e.g. Preibisch et al. 2003). Along 
these outflow cavities, FUV radiation may escape and 
irradiate the high density material at the border of the 
cavity. 

A surprisingly large amount of C0+ has been detected 
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in the high-mass star forming region AFGL 2591 by 
Stauber et al. (2007). Their radiative transfer calcu- 
lations indicate fractional abundances of order 10~^°, 
much higher than 10~^^ — 10"^'' predicted by dark 
cloud models. CO"*" has previously not been detected 
towards envelopes of young stellar objects (YSOs), but 
has been seen in photo dominated regions (PDRs, e.g. 
Jansen et al. 1995 or Fuente et al. 2003). Strong far 
UV (FUV) radiation in these regions heats the gas 
and drives a peculiar chemistry which enhances the 
abundance of CO"*". The detection of CO"*" in envelopes 
of YSOs is thus strong evidence for the feedback of 
the protostar on the envelope by high-energy irradia- 
tion. In this work, we will study CO"*" as a prototypical 
species enhanced by FUV irradiation in high- mass YSOs. 

Chemical models solve for the evolution of the 
abundances of molecules and atoms. They simulate 
a network of species reacting with each other. The 
network consists of different types of chemical reactions. 
In dark clouds for example, all FUV radiation is shielded 
by a large column density of dust and the chemistry 
is dominated by cosmic ray ionization and reactions 
between neutral and ionized species (e.g. Doty et al. 
2002). In the innermost part of a YSO envelope. X-rays 
may dominate the ionization (Stauber et al. 2005). 
They influence the chemistry mainly by secondary fast 
photoelectrons hitting molecules and atoms. X-ray 
induced chemistry is very similar to cosmic ray induced 
chemistry (Bruderer et al. 2009). In regions with 
a strong FUV irradiation, direct photoionization of 
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species with ionization potential < 13.6 cV drive the 
chemistry. The influence of FUV radiation on the en- 
velope of YSOs has been studied by Stauber et al. (2004). 

So far, a spherically symmetric structure of the 
envelope has been assumed for chemical models of 
envelopes of high-mass star forming regions (e.g. Viti 
& WiUiams 1999, Doty et al. 2002, Stauber et al. 
2004, 2005). While the abundances derived using this 
spherical model agree for many species with observations 
of specific sources, they fail to explain the amount of 
00+. Due to the attenuation by dust, FUV radiation 
cannot escape the innermost few hundred AU and 
the volume of gas where C0+ may be formed is too 
small in spherical models to reproduce the observations. 
Another possible formation mechanism of 00+ is X-ray 
irradiation. Unlike FUV photons. X-rays are mainly 
attenuated by geometrical dilution (cx r~^) due to their 
smaller absorption cross-section proportional to A^. 
X-ray radiation can thus penetrate much deeper into the 
envelope on scales of a few thousand AU. Protostellar 
X-ray luminosities of up to lO"*^ erg s^^ in the 1-100 
keV band arc observed (Prcibisch & Fcigelson 2005). 
For this luminosity however, the X-ray flux is too low to 
enhance C0+ sufficiently to match with observations. 

The interaction zone between outflow and infall has 
been proposed before as sites of anomalous chemistry 
and molecular excitation. Spaans et al. (1995) have 
explained the strong CO(6-5) emission in narrow lines 
observed towards many low-mass YSOs by heating of the 
cavity wall through protostellar FUV photons. Compact 
emission of HCO+ associated with the outflow wall has 
been detected by Hogerheijde et al. (1998) in the Class 
object L 1527. A possible explanation is mixing of 
partially ionized gas from the outflow with infalling 
material leading to a thin layer of dense and ionized 
gas, where species with a high dipole moment (e.g. 
HCO+ or HCN) can be excitet efficiently by collisions 
with electrons (Hogerheijde 1998). Detailed 3D radiative 
transfer studies by Rawlings et al. (2004) on the other 
hand prefer a scenario with an enhanced abundance of 
HCO+ due to shock liberation and photoprocessing of 
molecular material stored in ice mantles. Observational 
evidence for such an interaction layer in a high-mass 
star forming region are found by Codella et al. (2006) by 
studying the line profiles of different molecular species. 
Recently, van Kempen et al. (2009) flnd evidence for 
UV photons escaping through the outflow cones and 
impacting the walls for the low-mass protostar HH 46. 

In this work, we will introduce a chemical model that 
implements a non-spherical density structure taking an 
outflow cavity into account. A proper treatment of the 
FUV irradiated outflow walls requires two-dimensional 
chemical models with high spatial resolution. This 
is however computationally excessively expensive, 
especially if the temperature structure is calculated 
self-consistcntly with the chemical abundances. In 
the first paper of this series (Bruderer et al. 2009), 
we have introduced a new method for fast chemical 
modelling. A precalculated grid of chemical abundances 
depending on different physical parameters like the 
density, temperature or FUV flux is calculated and 



the abundances can then be obtained by interpolation 
in this database. Using different benchmark tests, we 
found very good agreement between fully calculated 
chemical abundances and interpolated values, while 
the interpolation approach is more than 5 orders of 
magnitudes faster. This speed-up allows to quickly 
construct detailed chemical models with a complex 
geometry. 

We construct a detailed two-dimensional chemical 
model of the high-mass star forming region AFGL 2591. 
The goal is to explain the high abundance of CO^ mea- 
sured by the JCMT single-dish telescope (Stauber et al. 
2007) as well as the non-detection by the Submillimeter 
Array (SMA) interferometer (Benz et al. 2007 and new 
data, reported in Sect. 4.4). In this paper, we assume 
a fixed geometry and consider only CO^. The next 
part of this series of papers is dedicated to a study 
of the influence of geometry for a larger sample of species. 

The paper is organized as follows: In the first section 
of the paper, we briefly discuss the chemistry of CO^. 
The next section introduces the two-dimensional model 
of AFGL 2591: We describe the modelling process, 
the radiative transfer of the FUV radiation and the 
calculation of the temperature structure. In section 
4, we present resiflts of the chemical model. Modeled 
fluxes and synthetic maps are compared to JCMT and 
SMA observations. 

2. CHEMISTRY OF CO+ 

The fractional abundance of C0+ depending on the 
temperature is given in Fig. 1 for different gas densities, 
X-ray fluxes [erg s~^ cm~^] and FUV irradiation in units 
of the interstellar radiation field (ISRF). A parcel of 
gas has been modeled for this figure using the chemical 
model described in the first paper of this series (Bruderer 
et al. 2009). The chemical model solves for the evolution 
of the abundances of a network of chemical species 
starting from initial abundances. We assume a chemical 
age of 5 X 10^ yr as suggested by Stauber et al. (2005). 
Chemical time-scales are very short in regions with 
strong FUV or X-ray irradiation, where a significant 
amount of C0+ is produced. Equilibrium conditions are 
thus reached very quickly and the chemical age does not 
influence the resulting abimdances. Photodissociation 
and ionization rates depend on the attenuating column 
density between the FUV/X-ray emission and the 
modeled parcel of gas. We assume an optical depth of 
T = 1 (FUV) and a hydrogen column density of 10^° 
cm~^ (X-rays) to derive the irradiating spectra. 

In gas with strong FUV irradiation, 00+ is effi- 
ciently produced by the reaction C+ -|- OH C0+ -|- 
H (Stauber et al. 2007, Fuente et al. 2006). Direct 
photoionization of carbon atoms accoimts for C^, while 
OH is produced by the photodissociation of water. At 
temperatures bc;low 100 K, H2O is frozen-out and the 
amount of OH reduced. Between 100 K and about 250 
K, the reaction O -|- H2 ^ OH -|- H does not proceed due 
to an activation barrier and OH is destroyed through 
photodissociation. The amount of OH increases with 
higher temperature and thus CO"*" is produced. At 
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very high FUV fields, C0+ is photodissociated and its 
abundance again decreases. 

In X-ray irradiated gas, ionization of CO through 
secondary electrons accounts for the production of C0+. 
The abundance of C0+ is thus constant with density 
and the fractional abundance ^ inversely proportional 
to the total density. In low density gas with high X-ray 
irradiation, the production of C0+ through C+ + OH 
can also be important as Figure 1 shows at a density of 
10^ cm~^ and an X-ray flux of 10 erg cm~^ s~^. In this 
regime, the fractional abundance of C0+ increases with 
temperature. 

CO"*" is destroyed by H2 leading to HOC"*" and 
IICO+ with equal branching ratios. In regions with very 
high FUV irradiation, H2 is photodissociated to atomic 
hydrogen and 00+ is destroyed by the reaction 00+ + 
R^CO + H+. 

To produce the observed fractional abundances of 
order 10"^" (Stiiubcr et al. 2007) by X-rays, a flux 
higher than 10 erg cm^^ s^^ is needed for a density of 
10^ cm~^. For an X-ray luminosity of 10"^^ erg s"""^ as 
found by Stauber et al. (2005) for AFGL 2591 and in 
the extreme case of no attenuation, this flux can only be 
achieved in the innermost 60 AU. In a spherical model 
of AFGL 2591, X-rays may enhance the abundance of 
C0+ up to a few times 10"", and Stauber et al. (2005) 
have thus suggested the molecule as a tracer for X-rays. 
X-rays are however not able to enhance the abundance 
to 10""'^'^ on larger distances to the source. The FUV 
flux of a young B star on the other hand is still sufficient 
at a distance of almost 10 000 AU to produce a fractional 
abundance exceeding 10~^° if no dust attenuates the 
FUV field. We thus consider C0+ to be a tracer for 
FUV irradiation rather than for X-rays. 

A short note concerning the combined influence of X- 
rays and FUV radiation: in the vicinity of a strong FUV 
source, CO is dissociated and X-rays cannot produce 
CO"'" efficiently. For example, the calculation of the pa- 
rameter study including X-rays at a density of 10^ cm~^ 
has been repeated in the presence of an FUV field of 10* 
ISRF. For this condition, the fractional abundance does 
not exceed 10~^* and we thus expect the fractional abun- 
dance in a low density and X-ray/FUV irradiated outflow 
region to be very low. 

3. AN AXI-SYMMETRIC MODEL OF AFGL 2591 

Figure 2 summarizes the process for an axi-symmetric 
model of AFGL 2591. We first assume a density struc- 
ture in Sect. 3.1. This density distribution is used in the 
following section to calculate the local FUV and X-ray 
field at every position of the model. The high-energy 
flux together with the density structure then enters the 
thermal balance calculation, which solves for the temper- 
ature structure. Since cooling and heating depend upon 
the composition of the gas, the temperature must thus 
be calculated self-consistently with the chemical abun- 
dances which are read out of the grid of chemical models 
presented by Bruderer et al. (2009). The abundance of 
C0+ is then obtained using the same approach. Syn- 
thetic maps are calculated using a raytracer and con- 



volved to the appropriate telescope beam for comparison 
with single-dish observations. To compare model results 
to interferometric observations, the synthetic maps are 
converted into visibility amplitudes and reduced in the 
same way as interferometric observations. 

3.1. Density structure 

One goal of this work is to study the difference in 
the modeled CO"*" line strength using a spherically 

symmetric ID and an axi-symmetric 2D model. For 
the 2D model we thus assume a density distribution 
following the power-law of the spherically symmetric 
ID model, except for an outflow- region (Fig. 3a). We 
implement n(H2) = no{ro/r) with no = 5.8 x 10^ cm~^ 
and To = 2.7 X 10"* AU, as proposed by van der Tak 
et al. (1999). 

For the separation between the outflow region and the 
envelope, we use the function 

^ " (l0 000tan2(a/2)) ' 

where z and r denote coordinates along the outflow axis 
and perpendicular to the outflow (in units of AU). The 
full opening angle a at 2; = 10 000 AU is assumed to be 
« 62°, approximately in agreement with the high resolu- 
tion mid-IR observations by Preibisch et al. (2003). The 
full opening angle of the outflow cone at z = 30 000 AU 
is about 40 degrees. The choice of a parabolic outflow 
cavity is supported by theoretical predictions (Canto 
et al. 2008) and observations (Velusamy & Langer 1998). 
It is also compatible with the mid-IR observations of 
AFGL 2591. 

Vibrationally excited CO at a velocity of —200 km 
s~^, emitted in outflow gas has been observed by van 
der Tak et al. (1999). Assuming an infall velocity in 
the envelope of order 10 km s~^, pressure equilibrium 
(pinWin = Pout f out) requircs a density ratio of 2.5 x 10~^ 
between the two regions. In the outflow region, the 
density of the power-law distribution is thus assumed 
to be reduced by this ratio. As in the spherical models 
of Doty et al. and Stauber et al., we do not model the 
innermost r:!250 AU. We make use of the symmetry and 
only model positive values of r and z. In this quadrant, 
the model consists of 300 cells along the outflow (2-axis) 
and 450 cells parallel to the midplane (r-axis) thus a 
total of 135 000 cells. The resolution is thus about 100 
AU X 67 AU. 



3.2. FUV and X-ray radiative transfer 

High-energy radiation with photon energy in the 
range of 6 - 13.6 eV (FUV radiation) and 0.1 - 100 keV 
(X-rays) penetrates molecular gas, ionizes molecules and 
heats the gas. The local FUV and X-ray field at every 
position are thus required to calculate the temperature 
structure and the chemical abundances. 

For the local X-ray intensity, we take into ac- 
count the distance r and the total hydrogen column 
density A^(Htot) [cm^^] between the X-ray source 
and the local position. For an X-ray luminosity Lx 
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Fig. 1. — Parameter study of CO+ : The fractional abundance is shown for a temperature range between 10 and 1500 K, different 
FUV fluxes between 0.1 and lO'^ ISRF at r = 1 in units of (top panels) and X-ray fluxes between 10"^" and 10 erg cm~^ for an 
attenuating hydrogen column density of 10^" cm~^ (bottom panels). Plots for hydrogen densities of lO'', 10® and 10^ cm~^ axe given. For 
compajrision, observed abundances are of order 10"^". 
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Fig. 2. — Flowchart of the modelling process for the axi- 
symmetric model of AFGL 2591. 

[erg s~i], the local intensity is thus proportional to 
Lx/47rr2 x cxp (-crphoto(i^) • iV(Htot)), where crphoto(i^) 
[cm^] is the photoionization cross-section at a photon 
energy E [eV]. We follow Stauber et al. (2005) and 
assume an X-ray emission by a thermal plasma at a 
temperature of 7 x 10^ K and a protostellar X-ray 
luminosity of lO''^ erg s^i. In their spherically sym- 
metric models, the attenuating column density in the 
innermost 250 AU is a free parameter found to be about 
3 X 10^^ cm~^. For the 2D model, the same value is 
adopted in the midplane (z = 0). Since density in the 
outflow is a factor of 400 lower, the column density 
between source and outflow walls can be of order 10^° 



cm~^, and photons with energies between 0.1 and 1 
keV are not absorbed. Unlike Stauber et al. (2005), 
we thus take into account these photons in the 2D model. 

The central FUV source is assumed to emit a black- 
body spectrum with a temperature of 3 x 10^ K and a lu- 
minosity of 2x 10*1 L0 (Stauber et al. 2004). At a distance 
of 200 AU, the FUV flux in absence of any attenuation 
corresponds to 2.3 x 10* times the ISRF. Unlike elastic 
compton scattering for X-rays (Stauber et al. 2005), scat- 
tering of FUV photons on dust grains can be important 
due to the similar wavelength of FUV radiation and the 
dust size. A Monte Carlo simulation in the vein of van 
Zadelhoff et al. (2003) calculates the local FUV strength 
(Appendix A). Photon packages starting at the central 
source are traced through the density model, while scat- 
tering and attenuation by dust is taken into account. For 
simplicity, we do not calculate the FUV spectra explic- 
itly, but deal with the scaling factor Gq in units of the 
ISRF. To read out chemical abundances from the grid of 
chemical models, the visual extinction r [Ay] is needed. 
We obtain the local extinction from 



(r) = - log 



(-^att)» 
(-^unatt)» 



(2) 



where /^^j is the local flux for a photon package i. 
The unattenuated flux /^natt '^^^Y geometrically 
diluted. Brackets denote averaging over all photon 
packages passing through a cell. We implement the 
cross-section on dust given in Draine (2003) for average 
Milky Way dust with Rv=3.1 and C/H=56 ppm in 
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Fig. 3. — Physical structure of the (axi-symmetric) 2D model of AFGL 2591: a.) density structure (Sect. 3.1). The border between the 
outflow and infall-region is indicated by the dotted line, b.) Attenuated FUV flux calculated by the FUV radiative transfer code (Sect. 
3.2). The gray/red arrow indicates the width of the outflow wall until Go is attenuated to 1 ISRF. c.) Temperature structure (Sect. 3.3) 
d.) Averaged value of (t) in units of Auv (Sect. 3.2). 



PAHs. This choice however has a minor effect on the 
model results, as photionization no longer dominates 
in regions, where a higher value of Ry « 5.5 is ex- 
pected. Since the adopted approach corresponds to 
the calculation at one FUV wavelength, we use the 
dust properties at a photon energy of 9.8 eV, in the 
middle of the 6 - 13.6 eV FUV band. For the conversion 
N(Rtot)/Av = 1.87 X 10^1 cm-2 (Bohlin et al. 1978) 
and an extinction cross-section of 1.29 x 10~^^ cm^ per 
H-atom at 9.8 eV, we obtain a conversion factor between 
the FUV and visual extinction of Auy/Ay = 2.4. 

For the density profile of Sect. 3.1, the calculated 
FUV flux (/att) and extinction (tuv) are given in Fig. 
3b and Fig. 3d. The attenuated flux shows the region 
where FUV radiation influences the chemistry, as this 
value enters the heating rates and the photodissociation 
rates of the form k — Gq ■ C ■ exp {—j ■ t^), with C 
and 7 fitting constants (e.g. van Dishoeck 1988) and 
Go the FUV flux in units of the ISRF. Note that the 
unattenuated flux (/unatt) has to be used together with 
the attenuation Tv to read out the abundances from the 
chemical grid. 



3.3. Temperature calculation 

Traditional PDR models (e.g. Tielens & HoUenbach 
1985) find the photoelectric effect of FUV photons on 
small dust grains and PAHs to be the dominant heating 
process in regions with strong FUV irradiation. In the 
dense outfiow walls of the 2D model, this process can 
easily heat the gas to temperatures above 250 K. For 
temperatures between 250 K and several hundred K, 
the fractional abundance of C0+ is a strong function 
of temperature (Sect. 2). It is thus important to carry 
out a detailed calculation of the temperature profile in 
the outflow wall. Input parameters are the gas density 
n [cm-^], the FUV held given by Gq [ISRF] and r [Ay] 
and the energy deposition per density by X-ray photons 
Hx [erg s"^]. 



To obtain the gas temperature Tgas, 



the thermal bal- 
ance has to be solved. We assume steady state condi- 
tions. Thus the equilibrium between the total heating 
(Ftot) and cooling (Atot) rate yields the temperature. 
Different physical processes contribute to the individual 
rates F^ and A^, where i denotes processes. The inner 
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energy per parcel of gas e [erg cm is given by 
Se 



St 



— Ttot ~ Atot 



= ^Tiin,GQ,T,Hx,xii^, . . .,N{CO), . . . , Tgas, Tdust) 

i 

- ^ A^{n, Go,T,Hx,Xii^,..., N{GO), Tgas, Tdust) ■ 

i 

Important examples for heating processes are H2 
collisional de-excitation following pumping through 
FUV photons and the photoelectric effect. The cool- 
ing rate consists of gas-grain collisional cooling and 
atomic/molecular line emission (e.g. the [CII] line at 
158 fim). An overview of the implemented cooling and 
heating processes is given in Appendix B. 

Heating and cooling rates depend on the compo- 
sition of the gas. For example the C+ atomic line 
cooling depends on the fractional abundance a;c+ of 
ionized carbon. To obtain a self-consistent solution, 
the composition of the gas is read out from the grid of 
chemical models. Because the temperature also enters 
as a parameter for the interpolation of the fractional 
abundances, the equation for the thermal balance has to 
be solved iteratively. 

Molecular or atomic line emission contributes to the 
cooling only if the photons can escape. A radiative 
transfer calculation is thus required to obtain proper 
cooling rates. For the atomic fine structure lines, we 
implement an escape probability code. As a simplifica- 
tion, the radiative transfer is not carried out in full 2D, 
but from the modeled point to the closest point of the 
outflow wall (Fig. 4). The probability (3 for a photon 
to escape from the outflow wall to the cavity is thus 
obtained from the molecular/atomic column density 
(e.g. ^(CO) for CO) along this line. 



outflow 
cavitv 




Fig. 4. — Schematic overview on the 2D modeL In the outflow 
waU, the simpHfied radiative transfer for heating and cooling is 
calculated between between the local position and the closest point 
of the outflow region. The column density (e.g. Af(CO) or A'(H2)) 
along this line is used to calculate the escape probability /3 of a 
photon. 



A fully self-consistent temperature model requires a 
two-dimensional radiative transfer calculation of the 
dust continuum radiation to obtain the dust temperature 
Tdust- As a simpliflcation, we use the dust temperature 
of the spherically symmetric model of Doty et al. 
(2002) for the region where the FUV radiative transfer 
calculation predicts an optical depth larger than 10. In 
regions with ry < 10, the dust temperature is obtained 
from the analytical expression given by Hollenbach 
et al. (1991). Using this approximation, the tempera- 
ture structure remains consistent with the spherically 
symmetric model for temperatures below 100 K. In the 
particular case of CO^, the relevant temperature range 



is T > 250 K, and the interpolation to the spheri- 
cally symmetric temperature profile is thus unimportant. 

4. RESULTS 

4.1. Directly irradiated outflow walls 

The adopted separation between envelope and outflow 
(z oc r^) yields a "flared" geometry that allows direct 
irradiation of the dense outflow walls by the central FUV 
source. This results in a much larger FUV irradiated 
volume compared to a "linear" separation (z oc r) where 
FUV photons can only penetrate the outflow walls if 
they are scattered on dust in the outflow (Fig. 5). 
Scattering is however relatively inefficient. 



Direct irradition 



Only scattered radiation 



,■>■ outflow wall 



envelope 



outflow 
-' Scattering .^.f 



outflow wall 



envelope 



Fig. 5. — Schematic view of the axi-symmetric model with conical 
and concave outflow cavity. 



What is the size scale of the outflow walls? We define 
the width of the outfiow wall by the region with Gq > 1 
along a cut of constant z. For example at z = 15 000 
AU, this corresponds to a thickness of about 1 130 AU. 
The outflow cavity thus allows long-range streaming 
of FUV radiation to large distances from the central 
source. Photon packages do not necessarily follow cuts of 
constant z. The straight ray path from the FUV source 
to the far end of the outflow wall at z = 15 000 AU 
propagates for 9 000 AU in the outflow wall. Scattering 
thus extends the outflow wall. For comparison, the 
Monte Carlo code was re-run with scattering switched 
off. This run yielded a width of the outflow wall of only 
about 470 AU at z = 15 000, about 2.5 times smaller 
than the width calculated including scattering effects. 
This demonstrates the importance of scattering for the 
strengths of the local FUV field. 

The width of the outfiow waU with Go > 1, 10 or 100 
ISRF depending on z is given in Fig. 6 along with the 
error due to the finite resolution of the model. An FUV- 
enhanced layer along the outflow with a width of a few 
hundered AU is thus produced by the protostellar FUV 
radiation. While the width defined by Gq = 1 increases 
with z, the layer with Gq > 100 has an approximately 
constant width of 280 ± 30 AU. This is a consequence of 
the adopted geometry: For the Gq > 100 layer, scatter- 
ing is less important than for the Go = 1 region. Unscat- 
tered radiation at the surface layer of the outfiow wall 
travels a longer way in the high density region for higher 
values of z due to the larger impact angle. This effect 
compensates with the fiux at the surface of the outflow 
wall oc l/d^, where d is the distance from the FUV source 
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to the surface of the outflow wall, and the attenuation 
decreasing with l/d for the implemented density profile. 



^ 1250 
I 1000 



T = 100 K 



ji?'* Go = 1 ISRF 



x5l 



,4}'" .m5"5*"'"'"""Go = 10 ISRF 
Go = 100 ISRF 



5 • lO' 1 • 10* 1.5 • 10* 2 • 10* 2.5 • 10* 

z-Axig [AU] 

Fig. 6. — Width of the outflow waU depending on the distance 
along the outflow axis. The figure sliows the width defined by the 
region with temperature > 100 K (line) and with FUV flux larger 
than 1, 10 or 100 times the ISRF (error bars). 

Only a small fraction of the envelope in a thin layer 
on the outflow wall is heated by FUV photons to 
temperatures above 100 K. The temperature along cuts 
of fixed values of z is shown in Fig. 7. For z < 5000 
AU, the temperature at the surface exceeds 1000 K, 
consistent with the surface temperatures obtained from 
FDR models for similar conditions (e.g. Kaufman et al. 
1999). 

The width of the region with T > 100 K in the 
outflow wall is approximately constant with z (Fig. 6) 
as the Go = 100 layer. In the region with T » 100 
K, the heating is dominated by the photoelectric 
effect TpE oc n ■ Go, with n being the density. The 
main coolant is atomic oxygen which is approximately 
thermalized and the cooling rate thus Aq oc n. As 
found in the calculation, the gas temperature is thus 
approximately independent of density and thus follows 
Go. 

An overview of the mass and volume of different 

regions of the axi-symmetric 2D model compared to 
the spherical ID model is given in Table 1. Indeed, 
the volume heated to T > 100 K is more than two 
orders of magnitude larger in the 2D model, while 
the mass of the same region is only « 8 times larger 
due to the adopted density gradient. Similarly, the 
volume irradiated by Gg > 1 ISRF is almost 4 orders 
of magnitude larger in the 2D model, while the mass 
of the gas is about 2 orders of magnitude larger. 
We note that size of the region with Go > 1 ISRF 
in the spherical ID model heavily depends on FUV 
attenuating column density in the innermost 250 AU. 
For Table 1 the extreme case of no attenuation be- 
tween the protostar and a radius of 250 AU was assumed. 

Physical parameters like the cavity shape or the 
density ratio between infalling envelope and outflowing 
gas enter the modelling. A complete study is beyond the 
scope of this study and will be presented in the third 
paper of this series (Bruderer et al., in prep.). Here we 
discuss some general trends. A large surface directly 
irradiated by protostellar FUV radiation is required 
for a significant volume and mass of the dense and 



FUV heated material in the outflow wall. By changing 
the outflow angle a by ±20 degrees or the separation 
between outflow region and envelope to z oc r" with 
a = 1.5 — 2.5, the volume and mass with T > 100 K 
vary only by a factor of about two. On the other 
hand, a higher density in the outflow region affects the 
model much more due to the higher extinction of FUV 
radiation. Increasing the density ratio nout/^in from 
2.5 X 10-3 to 10-2 reduces the mass of the T > 100 K 
region by a factor of 4, while the volume decreases by 
about an order of magnitude. These results for the 
relatively insensitivity with respect to angle and shape, 
and relative sensitivity to density, underline one of the 
key results of this paper namely the importance of 
direct irradiation on the chemical composition of outflow 
cavity walls. 

The calculation of the FUV flux and gas temperature 
in the axi-symmetric 2D model suggests a region with 
temperature T < 100 K but FUV irradiation Go > 1 
of potentially FUV irradiated ices. While a detailed 
discussion of these ices is beyond the scope of this 
paper, it is interesting to note two things. First, such 
a region exists due to the direct irradiation allowed by 
the concave outflow walls. Second, the potential mass in 
such a region is 1.1 Mq. Together, these facts suggest 
that further exploration of irradiated ices is warranted. 



4.2. CO^ abundance 

Given the physical parameters for an axi-symmetric 
2D model of AFGL 2591 defined previously, the grid of 
chemical models (Bruderer et al. 2009) can now be used 
to quickly read out the abundances of C0+. Since no 
outflow features have been observed in lines of C0+ by 
Stauber et al. (2007), we assume there is no C0+ in the 
outflow cavity. Also, the discussion in Sect. 2 indicates 
a very low fractional abundance of CO"*" in the outflow. 

The fractional abundance of CO"*" is given in Fig. 7 
along with the temperature and the density for cuts 

parallel to the midplane of the model. In a thin layer at 
the outflow wall, CO"^ is enhanced. This layer matches 
with the FUV heated gas at temperatures above 250 K. 
The C0+ layer in the midplane {z = 0) is considerably 
larger due to the scattering of FUV radiation from both 
outflows along the positive and negative z-axis. At the 
surface of the outflow wall, the fractional abundance of 
CO"'' exceeds 10^^° for values of z up to a few thousand 
AU. At z = 25 000, the fractional abundance of C0+ is 
still a few times IQ-^^. 

Since the chemical time-scales of C0+ in regions with 
abundance larger than a few times 10^^^ are shorter 
than 100 yrs, the evolutionary age of the source does 
only affect the modelling results if the FUV production 
of the protostar has not set in yet or the outflow cavity 
is not large enough yet. While direct evidence for pro- 
tostellar FUV radiation of the sourc;e c;annot be given, 
the mid-IR observations of the source by Preibisch et al. 
(2003) suggest an outflow cavity with size of at least 
10 000 AU along the outflow axis. For the measured fiow 
velocity of a few hundered km s^^ (van der Tak et al. 
1999), the expansion to 10 000 AU would require < 500 
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TABLE 1 

Volume and mass of different regions of the axi-symmetric model compared to the spherically symmetric model. 



Region 


Mass 


Vol 


ume 




[Mq] 


[AU^] 


Axi-symmetric 2D model: 






1013 


Outflow 


0.014 


1.1 


Envelope 


44.7 


1.0 


1014 


Envelope (T > 100 K) 


0.53 


8.2 


lOii 


Envelope (Go > 1) 


1.65 


3.0 


1012 


Envelope (Go > 1 and T < 100 K) 


1.1 


2.2 


1012 


Spherical ID model: 








Envelope 


50.3 


1.1 


1014 


Envelope (T > 100 K) 


0.07 


6.5 


■ 109 


Envelope (Go > 1) 


0.01 


4.8 


■ 108 



yrs much less than the infered age since YSO formation 
estimated to be a few times 10^ yrs (Doty et al. 2002, 
2006). 

For comparison, results of the spherical ID model 
are given in the top left of Fig. 7. In this model, 
FUV radiation of the central source cannot escape the 
innermost part and the abundance of CO"'" is thus very 
small at fractional abundances below 10~^^. Outflow 
walls in the 2D model increase the surface irradiated by 
FUV and thus enhance the total amount of CO"*". 

Mixing between warm and ionized outflow material 
with the envelope can also induce a particular chemistry 
in a thin layer along the outflow. Charge-exchange re- 
actions in this mixing layer enhance the abundance of 
CO"*". The width of the mixing layer is however small 
due to the high electron abundance leading to short re- 
combination time-scales. While a detailed study (Taylor 
& Raga 1995, Lim et al. 1999 or Nguyen et al. 2000) is 
beyond the scope of this work, a simple model is given in 
Appendix C. It predicts the amount of CO""" produced in 
this mixing layer to be enhanced, however insuSiciently 
to explain observations. 

4.3. Comparison with JCMT observations 

Two rotational lines of CO"*" with Av 4 km s""^ have 
been detected by Stauber et al. (2007) at the center 
position of AFGL 2591 using the JCMT. We follow 
their arguments and assume the lines to be detected. 
The observed line fluxes are given in Table 2 along with 
molecular data of CO"*" taken from the JPL database 
(Pickett et al. 1998). The Einstein- A coefficients are 
scaled to a dipole moment of 2.77 D (Rosmus & Werner 
1982). 

Synthetic spectra arc calculated to compare the 
abundance modelling results to these observations. The 
excitation of CO"*" is needed to model the line flux. 
However, the mechanism is unclear since CO"*" is more 
more likely to be destroyed than excited in collision 
with II2. However, electrons and atomic hydrogen arc 
abundant at the surface of PDRs and can excite C0+ as 
discussed in Andersson et al. (2008). Possibly, C0+ is 
formed in an excited state. In this work, we follow 
Stauber et al. (2007) and assume a fixed excitation 
temperature Tex- The effects of collisional excitation by 
electrons and atomic hydrogen and excited formation 
are being studied by Stauber & Bruderer (2009). For 



the transitions discussed in this section and assuming 
an excited formation, they find a fiux increase by a 
factor of only w 3 compared to the peak LTE flux at an 
excitation temperature of w 30 K. 

To constrain rotational temperatures, transitions for 
different upper levels are necessary. However, only 
tentative detections of C0+ in YSOs for Jup > 14 are 
reported (Ceccarelli et al. 1998). Wc will thus follow 
Black (1998) and assume Tox=30 K, slightly higher than 
Stauber et al. (2007) (20 K) or Fucnte et al. (2006) (10 
K). For an upper energy level Eu, the level population is 
proportional to e-^''/^'^«/Q(Tex), with Q{T) being the 
partition function. Using Tox = 30 K, the population 
of the J = 3 levels is approximately maximized. The 
emission of optically thin lines scales with the population 
of the upper level. Compared to Tex = 30 K, the line 
flux is thus about about 20% weaker for Tex = 20 K and 
a factor of 3.5 for Tex = 10 K or T^^ = 300 K. An exci- 
tation temperature of a few hundred K is expected for 
collisional excitation by electrons and atomic hydrogen 
(Andersson et al. 2008). We conclude that a different 
excitation temperature, within the expected range or 
assuming an excited formation of CO"*", changes the 
modeled line flux less than an order of magnitude. 

Synthetic line fluxes are obtained from the solution of 
the radiative transfer equation. The two lines considered 
in this work have a very low optical depth with t < 0.08 
for rox=30 K at an assumed line width and molecular 
column density of 4 km s^^ and 9 x 10^^ cm~^, respec- 
tively. This column density corresponds to an emitting 
region with a diameter of 60 000 AU at a density of 10 
cm~^ with a CO"*" fractional abundance of 10~^°. For our 
application, we can thus neglect self absorption and the 
integrated flux is obtained from the velocity integrated 
radiative transfer equation 

with V being the line frequency, Qu the statistical weight 
of the upper level, Aui the Einstein-A coefficient [s^^] 
and n{s) the C0+ density [cm^'^]. The integral is 
taken along the line of sight (LOS). To have comparable 
results to older literature, we assume a distance of 1 
kpc to AFGL 2591, following Benz et al. (2007) and van 
der Tak et al. (1999). A larger distance of 1.7 kpc as 
suggested by Schneider et al. (2006) would however not 
change the main conclusions of this work. The line flux 
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Fig. 7. — CO"*" in the two-dimensional axi-symmetric model of AFGL 2591. Slices at different values of z are shown. The x-axis 
corresponds to the distance along the cut to the outflow region. The plots give the fractional abundance a:(CO"*") = n(CO+)/n(H2), the 
temperature and the density in solid, dashed and dotted lines, respectively. The figure in the top-left shows a spherical model of the same 
source. The gray rectangle indicates the innermost 250 AU where no modelling has been performed. 

of the JCMT. The effect of different inclination angles 
is clearly visible in unconvolved maps, but the con- 
volved maps are almost independent of inclination angle. 



of the spherical ID model calculated using this method 
agrees within a few percent with results from SKY of the 
RATRAN radiative transfer code (Hogerheijde & van 
der Tak 2000). 

Maps of the integrated line flux from the 2D model 
at different inclination angles are shown in Fig. 8. 
The angular resolution of the images is 0.05". We give 
only results for the 3| - 2| transition at 354.014 GHz 
in this section. The modeled flux of the other fine 
structure line can be obtained by scaling with 0.7, the 
ratio of the Einstein-A coefficients and the statistical 
weights. The modeled maps are convolved with a 14" 
FWHM Gaussian to account for the angular resolution 



The predicted line fluxes for the JCMT at the center 
position are given in Table 3. They only vary between 
0.16 and 0.14 K km s"^ for an inchnation of 0° and 90°, 
respectively. This reflects a slightly larger fraction of the 
outflow walls in the telescope beam for low inclination 
angles. The difference is however small since most 
emission is from regions close to the center. Indeed, 
about between 15% and 30% of the modeled JCMT 
flux stems from within a radius of 1" from the center, 
corresponding to 1 000 AU at the adopted distance. An 
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TABLE 2 

Molecular data of CO+ and integrated line flux at the center position of AFGL 2591 observed by the JCMT (Stauber 

ET AL. 2007). 



Transition 


Frequency 
[GHz] 


Aul 

[10-3 s-1] 


[K] 


Line Flux 
[K km s-1] 


3| - 2| 

'il _ 9^ 
■^2 ^2 


353.741 


1.58 


33.94 6 


0.5 


354.014 


1.7 


33.99 8 


0.27 



inclination angle of 30° was suggested by van der Tak 
et al. (1999) and is thus displayed in this and the next 
section. 

The line flux of the C0+ hne at 354.014 GHz as ob- 
tained from the spherically symmetric ID model is with 
1.1 X 10~^ K km s~^ about 3 orders of magnitude weaker 
than observed. On the other hand, the results of the 
axi-symmetric 2D model agree to within a factor of two 
which wc consider as a good agreement given all un- 
certainties entering the modelling. The 3| - 2| line 
at 353.741 GHz however disagrees by about a factor of 
5, as the observed line ratio of 1.85 (=[353.741 GHz] 
/ [354.014 GHz]) is larger than the ratio of 0.7 as pre- 
dicted by the models. Possible explanations arc a dif- 
ferent excitation mechanism of the fine structure levels 
or line blending of the 353.741 GHz line, e.g. by the 
33S02(JKpX„ = 194,16 ^ 193,17) line at 353.741 GHz 
(Stauber et al. 2007). To explain the line ratio by a 
different excitation iiiechanism, a grossly different exci- 
tation temperature would however be required, as dis- 
cussed earlier in this section. 

4.4. Comparison with SMA observations 

Interferometric observations of the two C0+ lines in 
AFGL 2591 have been carried out using the Submil- 
limeter Array (SMA)^. Benz et al. (2007) have used 
the extended configuration of the array with projected 
baselines covering the range between 38 - 214 kA (32.2 
m - 181.2 m). New observations have been carried out in 
the compact configuration on April 14 2006. These data 
cover a projected baseline range of 10.6 - 82 kA (9 m - 
69.4 m). The frequency setting is the same as for the 
extended array observation, covering the range between 
342.6-344.6 GHz in the LSB and 352.6-354.6 GHz in 
the USB. Data of this observation will be discussed in a 
paper in preparation (Bruderer et al. 2009). 

The RIVIS noise per velocity bin of 0.5 km is 0.46 
Jy Beam"^ (extended array), 0.21 Jy Beam-^ (compact 
array) and 0.18 Jy Beam"-*^ when data of both arrays are 
combined using natural weighting. The better weather 
conditions and shorter baselines result in a higher 
sensitivity of the compact array observations. The 
detection limit for the integrated line can be estimated 
from Icr = I.2VSY ■ 6vT,^^, where T^ms is the RMS 
noise per frequency bin, (5V ~ 4 km s~^ the expected 
linewidth and 5v = 0.5 km s~^ the velocity resolution 
(Maret et al. 2004). The factor 1.2 accounts for the 
uncertainty of the absolute calibration of about 20%. 

® The Subinillimeter Array is a joint project between the Smith- 
sonian Astrophysical Observatory and the Academia Siniea Insti- 
tute of Astronomy and Astrophysics and is funded by the Smith- 
sonian Institution and the Academia Sinica. 



For a 3(7 detection, an integrated flux of 2.3 Jy km s~^ 
(extended array) and 0.91 Jy Beam"^ km s~^ (compact 
and extended arrays combined) is thus necessary. 

The C0+ hne at 354.014 GHz has not been de- 
tected at a 3(7 level in either configurations of the 
SMA. The 353.741 GHz line is not detected in the 
extended configuration observations. In the data of 
the combined array however, a line is found at ap- 
proximately 353.741 GHz. Is the C0+ line blended by 
^^S02{JKpKo = 194,16 193,17) at the same frequency? 
Several other lines of ^^S02 are within the observed 
frequency setting, but are not detected despite their 
presumably higher intensity due to a larger Einstein-A 
coefficient, lower critical density and upper level energy. 
G0+ is thus more likely blended by an unidentified line. 
The blended line has a velocity integrated line flux of 4.9 
Jy Bcam-'^ km corresponding to a 25a detection. 
For the predicted hne ratio ([353.741 GHz] / [354.014 
GHz]) of 0.7 and the ratio of 1.85 as observed by the 
JCMT, the C0+ hne at 354.014 GHz should have been 
detected. We thus conclude C0+ is not detected by the 
SMA observations. 

Is this non-detection consistent with the modeled line 
flux of the 2D model? To answer that question, synthetic 
maps (Fig. 8) are converted to visibility amplitudes for 
the (u, v) coverage of the SMA observation using the 
MIRIAD^ task uvmodel. The simulated visibilities are 
then reduced in the same way as SMA observations: 
they are inverted and the clean algorithm is applied 
on the resulting maps. Simulated velocity integrated 
maps for inclination angles of 30° and 90° are presented 
in Fig. 9. The maximum flux of the simulated maps 
for the extended array is 0.17 Jy Beam^^ km s^^ at 
an inclination of 70°. For the (u,v) coverage of the 
combined array, the peak flux is 0.47 Jy Beam"^ km s~^ 
at an inclination angle of 60°. These integrated fluxes 
correspond to less than Icr and the non-detection is 
thus consistent with the models. As the simulated SMA 
maps in Fig. 9 show, this non-detection is explained by 
the limited sensitivity of our observations rather than 
by over-resolving effects of the interferometer. 



5. CONCLUSIONS AND OUTLOOK 

We have used the grid of chemical models introduced 
by Bruderer et al. (2009) to construct a detailed two- 
dimensional model of the high-mass star forming region 
AFGL 2591. The spherically symmetric ID models by 

Doty et al. (2002) and Stauber et al. (2004, 2005) have 
been extended with a low-density outflow region allowing 

http;//sma-www.cfa.harvard.edu/miriadWWW/manuals/ 
manuals.html 
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Fig. 8. — Synthetic maps of the velocity integrated CO+ hne at 354.014 GHz (3^ — 2^) for different inchnation angles of the source. 
The inclination angle is shown with the gray/red arrow onto a symbolic density map in the bottom-left corner of each map. Two plots in 
the lower right corner give maps convolved the a 14" JCMT beam. 

region. In the 2D model, FUV irradiates a larger surface 
and thus a larger volume compared to spherically sym- 
metric models. The model is used to simulate the line 
fluxes of C0+ as a prototypical FUV tracer. The main 
conclusions of this work are: 

1. The existence of concave outflow cavity walls al- 
lows the efficient, long-range streaming of FUV ra- 
diation to large distances from the central source 
(Sect. 4.1). 

2. A thin FUV-enhanced layer in the outflow is pro- 
duced, having a thickness of a few hundreds of AU, 
and Go > 100 ISRF (Sect. 4.1). 

3. A large mass of FUV irradiated material - more 
than a solar mass - can reside in the outflow walls. 
This combined with the large extent can lead to 
a situation where the molecular emission from the 
outflow walls dominates the total line flux in single- 
dish observations (Sect. 4.1, Sect. 4.3) 

4. Detailed 2D modelling of C0+ including these out- 
flow walls is consistent with single-dish (Sect. 4.3) 
and interferometric observations (Sect. 4.4), while 
the flux in a ID model is orders of magnitude too 
low. This indicates the need for multidimensional 
chemical models for the interpretation of FUV sen- 
sitive molecules. 
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Fig. 9. — Simulated SMA observations [Jy Beam ^ km s ^] 
(gray/red contour lines) overlaid on modeled velocity integrated 
maps [K km s~^] (black/white map). Negative fluxes are shown by 
dashed contour lines. The half-power beam is given in the bottom 
left corner. For comparison, a 3cr detection with the SMA requires a 
flux of 2.3 (only extended array) / 0.91 (compact -|-extended array) 
Jy Beam^^ km s^^. 

protostellar FUV radiation to escape from the innermost 
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TABLE 3 

Modeled and observed maximum velocity integrated fluxes of the CO+ 3| — 2| line at 354.014 GHz. 



Inclination 


SMA (Ext.'') 


SMA (Both.'") 


JCMT 


[°] 


[Jy Beam" ^ km s" ^] 


[Jy Beam^^ km s~^] 


[K km s-i] 


Axi-symmetric 2D model: 









0.05 


0.23 


0.16 


10 


0.06 


0.25 


0.16 


20 


0.05 


0.26 


0.16 


30 


0.06 


0.30 


0.15 


40 


0.06 


0.38 


0.15 


50 


0.11 


0.44 


0.15 


60 


0.16 


0.47 


0.14 


70 


0.17 


0.45 


0.14 


80 


0.16 


0.44 


0.14 


90 


0.16 


0.43 


0.14 


Spherical ID model: 








< 0.01 


< 0.01 


1.1 • 10-4 


Observed: 


< 2.3=^ 


< 0.91^= 


0.27" 



''extended configuration 
'"compact configuration 
"^Not detected: 3a upper limit 
''Detected 

5. The strong gradients in density and FUV flux re- 
quire an accurate calculation of the gas tempera- 
ture structure and inclusion of scattering of FUV 
radiation in the outflow walls, which can extend 
the width of the FUV enhanced region significantly 
(Sect. 4.1). 

6. A region of high FUV (Go > 1 ISRF) and low tem- 
perature {T < 100 K) leads to the strong possibil- 
ity of ice mantle processing by FUV photons in the 
outflow walls (Sect. 4.1). 

7. The abundance of C0+ is also enhanced in a mix- 
ing layer between the ionized and warm outflow 
and the envelope. Mixing alone can however not 
explain the observed amount of C0+ (Sect. 4.2, 
Appendix C). 



We thank Kaspar Arzner, Pascal Stiiubcr and Jes 
J0rgensen for useful discussions and an anonymus referee 
for his/her valuable comments. Michiel Hogerheijde and 
Floris van der Tak are acknowledged for use of their 
RATRAN code. This work was partially supported 
under grants from The Research Corporation and the 
NASA grant NNX08AH28G (SDD). Astrochemistry in 
Leiden is supported by the Netherlands Research School 
for Astronomy (NOVA) and by a Spinoza grant from 
the Netherlands Organization for Scientific Research 
(NWO). The submillimeter work at ETH is supported 
by the Swiss National Science Foundation grant 200020- 
113556. 

Facilities: SMA, JCMT 



This paper shows the application of the grid of chemi- 
cal models (Bruderer et al. 2009) for the construction of 
a multidimensional chemical model of a YSO envelope. 
This interpolation method simplifies the construction of 
two or three dimensional chemical models considerably 
and the gain in speed allows the self-consistent calcu- 
lation of the temperature structure with the chemical 
abundances (Sect. 3.3). Our work shows that such 
detailed modelling including the; influence of high-energy 
irradiation and the geometrical shape of the envelope 
can be necessary even to explain spatially unresolved 
single-dish observations. 

Some hydrides (e.g. SH+ or NH+) will be observable 
for the first time using the upcoming HIFI spectrome- 
ter onboard Herschel. They are expected to be chemi- 
cal tracers of FUV/X-ray radiation (Stauber et al. 2004, 
2005, 2007) and further study will be carried out to in- 
vestigate the influence of the geometry on these species 
for the interpretation of the Herschel data. Direct imag- 
ing of the outflow region will be available with the At- 
acama Large Millimeter/submillimeter Array (ALMA). 
These high spatial resolution observations will deliver ad- 
ditional constraints on the geometry. 
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APPENDIX 

A. A MONTE CARLO CODE FOR FUV RADIATIVE TRANSFER 

In the first appendix, we discuss the implementation of the radiative transfer calculation of the FUV radiation. 
Unlike the work of van Zadelhoff et al. (2003), our calculation is implemented on a three dimensional Cartesian grid 
and thus already able the solve three dimensional problems. To improve convergence, points with constant z and 
r = -\- "ip- are averaged for the axi-symmetric model in this work. To save memory, only cells with positve values 
of X, y and z are stored, but photon packages can of course still travel from one quadrant to another. The flow of the 
calculation is 

1. A photon package is started from the FUV source at an arbitrary direction with a luminosity of ^package = 
isourcc/(8 • -/Vphoion) [crg s~^], whcrc Lsouicc IS thc luminosity of the FUV source in the 6 - 13.6 cV band. A^phton 
is the number of photon packages in the simulation run. The factor 8 takes into account that only cells with 
positive values of x, y and z are stored. 

2. As the photon package travels along a straight line, ^package is attenuated by exp (— rioc) with rioc = 
(Text n{x, y, z) As, the local FUV optical depth. The optical depth depends on (Text, the extinction cross-section, 
n{x, y, z) the local density and As the step size of the code. The total extinction, that a photon package suffers 
from thc source to thc current position rtot = ''"loc is stored. To avoid discretization errors, especially close to 

thc source, wc use thc analytical expression for thc local density. 

3. If the photon package passes the border of a cell, the local intensity is updated. The attenuated intensity is 
obtained from 

J ( s _ \ ^ J [nphoton • Ocelli 

-'att V^) 2/, j — / J -^package ^ i 

where the sum is taken over all photon packages passing the cell. The surface of the cell wall which the photon 
package just passed is given by A [cm^] and the directions of the photon package and the normal of the cell 
separation is given by fiphoton and nceii, respectively. 

4. The unattenuated intensity /unatt is calculated using the same expression as 7att , except I/package is replaced 
by ^package X exp (+rtot)- 

5. If a random number [0,1] is larger than exp (— Tgcat): with Tgcat = CTcxt ri{x, y, z) As, thc scattering-optical depth, 
a new direction of the photon package is chosen following the implemented phase function (/'(A). For our work, 
we use the tabulated function of Draine (2003) for an average Milky Way dust with Ry = 3.1 and C/H=56 
ppm in PAHs. 

6. For each photon package, steps 2-5 are repeated, until ^package is smaller then a certain threshold or the package 
left the simulated area. 

7. The code propagates 10"'' photon packages and then averages for the attenuated intensity /att in the axi-symmetric 
model. Three subsequent solutions are used to calculate the signal-to-noise ratio (SNR), defined by 



SNR"^ = max 



^att (-^a 



Only cells with an averaged intensity larger than lO^"' ISRF are considered. Convergence is reached, if the SNR 
exceeds 50 in all cells. This means less than 2% difference to the average. The application in this work requires 
9 X 10^ photon packages to achieve this accuracy. 

As a benchmark test, the code is re-run for the density distribution in Sect. 3.1, but with scattering switched off. 
In this way, an analytical solution can be compared to the result of the code. For the analytical solution, only the 
distance r and column density NiYi) [cm~^] to the FUV source are needed to derive the local FUV fiux [erg s~^ cm~^] 

^^°^ = ^xexp(-ruv) 

with the FUV luminosity Lyv = 4.2 x 10^"^ erg s^^ and the FUV optical depth tuv = 2.4(7V(H) [cm~^])/1.87 x lO^i. 
Thc conversion factor of 1 ISRF = 1.6 x 10^^ erg s^^ cm^^ yields thc flux Go in units [ISRF]. Thc agreement between 
the code and the analytical calculation is tested for points with Go > 10~^. A very good agreement is found, with 
deviations less than 50%. 



B. THERMAL BALANCE CALCULATION 

In this appendix, references for the heating and cooling rates used in Sect. 3.3 are given and the temperature 
balance calculation is benchmarked with FDR models. 
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HEATING RATES 

Photoelectric heating: The heating rate of photoelectrons from small dust grains and PAHs is Fp^; = 10~^^ e n Go erg 
cm~^ s~^, with the density n, the attenuated FUV field Gq and an efficiency e. The efficiency depends on the 
ionization fraction of the grains, since the work function increases for ionized dust grains. Bakes & Tielens (1994) 
give an analytical fit for e as a function of the ratio 7 between ionization and recombination rate (7 = Go^/T/ne, with 
the electron density Ug and the temperature T). 

H2 heating: Different heating processes involving molecular hydrogen arc considered: (i) CoUisional dc-cxcitation of 
vibrationally pumped H2 by FUV photons (FUV-pumping). (ii) Formation heating: H2 forming on dust releases 
part of the binding energy to the gas. (iii) Photodissociation of H2 heats the gas. The rates for process (i) and (iii) 
depend on the local FUV field. Since line absorption is responsible for the pumping and dissociation of the molecule, 
we reduce the local FUV field by the self shielding factor /3(t) (Draine & Bertoldi 1996). Rates for (i) and (ii) are 
implemented by the analytical expression of RoUig et al. (2006). The rate of process (iii) is taken from Meijerink & 
Spaans (2005). 

C ionization: Photoionization of atomic carbon, C + 7uv ^ C+ + e~, releases an energy of 1.06 eV to the gas and 
thus contributes to the heating. We follow Tielens (2005) for the implementation. 

X-ray heating: Fast photoelectrons produced by X-ray photons lose part of their energy through Coulomb interaction 
with thermal electrons. The heating rate is given by Fx = r]nHx, where Hx [erg s~^] is the energy deposition of 
X-rays in the gas (c.f. Stauber et al. 2005) per density. The efficiency factor r] depends on the H/H2 ratio and the 
ionization fraction. An analytical fit for r] to detailed calculations is given in Dalgarno et al. (1999). 

Cosmic ray heating: For low degrees of ionization x < 10~^, about 9 eV per primary ionization through a cosmic ray 
particle is used to heat the gas. The heating rate is thus Fcr = 1.5 x 10~^^ erg cm~^ s~^. 

COOLING RATES 

Atomic fine structure lines: Forbidden fine structure lines of O, C and are important coolants at the surface of 
the FUV irradiated zone. We consider the [01] 63 fim, [01] 146 fim, [CI] 369 /xm, [CI] 609 /xm and [CII] 158 /zm hues. 
The cooling rate is obtained from 

Aline = ^ul hVul 13 {TuI ) n„ (X) , 

with the Einstein-A coefficient A^i [s^^] of the transition u ^ I, v^i [s^^] the line frequency, t„; the optical depth, 
(Hjui) the escape probability and n„(X) [cm~'^], the density of the species X in the excited level u. An escape 
probability formalism (e.g. Tielens 2005) is used to calculate the level population iteratively. CoUisional excitation by 
H, H2 and electrons is taken into account - the abundance of the collision partners is read out of the chemical grid. 
We use the same collision rate coefficients as Meijerink (2006). The optical depth r is calculated using the molecular 
column density to the outflow cone (Sect. 3.3). 

H2 line cooling: Vibrational lines of H2 can contribute to the cooling of the gas. Due to the large gap between the 
ground state and the first excited state, corresponding to about 6 000 K, we use the two level approximation given in 
RoUig et al. (2006) as a simplification. Again, self shielding is taken into account for the radiative pumping by FUV 
photons. 

Gas-grain cooling /heating: The difference in temperature between Tgas and the Tdust leads to a transfer of heat. This 
can be cooling close to the FUV irradiated zone, where Tjust < Tgas, or heating deeper in the envelope. The rates are 
proportional to T^ust — ?gas- We implement the results of Hollenbach & McKee (1989) for a minimal grain size of 10 

A. 

Cooling by CO and H2O: Molecules can contribute to the gas cooling by rotational lines at low temperature and 

ro-vibrational lines at higher temperature. We include line cooling by CO and II2O using the fitted rates of Neufeld 
& Kaufman (1993) (for T > 100 K) and Ncufcld ct al. (1995) (for T < 100 K). The rates depend on the molecular 
column density, the temperature and the density of the collision partner. Their work considered excitation by H2. 
Excitation by electrons and atomic hydrogen is taken into account by the approximation given in Yan (1997) and 
Meijerink & Spaans (2005). The molecular column density is obtained in the same way as for the cooling through 
atomic fine structure lines. Cooling by isotopes (-"^"^CO, C^^O, 112^0) can be important due to the smaller optical 
depth. Isotope ratios by Wilson & Rood (1994) are used to scale the molecular density and the column density. 

Recombination: Recombination of electrons with grains and PAHs can cool the gas. We implement the fit-results 
by Bakes & Tielens (1994) for our calculation. This rate also depends on the ratio 7 between ionization and 
recombination rate, due to increased Coulomb interaction in gas with a high ionization fraction. 

Lya emission: Atomic hydrogen excited by electron collisions emit Lya photons and thus contribute to the cooling. 
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This process is only eflicient at temperatures higher than a few thousand K. We use the cooling rate given in Sternberg 
& Dalgarno (1989). 

BENCHMARK TEST 

We compare our calculation of Tgas to PDR codes, in a similar situation as model V4 in the PDR comparison study 
by Rollig et al. (2007). The test consists of a plane parallel slab at a density of 10^-^ cm~'^ and a FUV irradiation 
of X = Go/1.71 = 10^. We assume an X-ray flux of about 0.04 erg s^^ cm^^ consistent to the AFGL 2591 model at 
z = 1000 AU assuming an X-ray luminosity of Lx — lO'^^ erg s~^. Differences in the resulting gas temperature, 
when X-rays are switched off, are however negligible. 



J a.} Gafi and Dust Temperature b.) Fractiorral Abundances 




Alteimation r [Av] Attemiation r [Av] 

Fig. B1. — Benchmark test of the heating/cooling calculation, a.) Result of the thermal balance calculation described in Sect. B. The 
gray shaded region gives the range of the PE)R comparision study by Rollig et al. (2007) (one standard deviation, mean=average=dotted 
line), b.) Fractional abundances of species relevant to the thermal balance calculation, c.) Heating rates, the line on top gives the total 
rate, d.) Cooling rates, the line on top gives the total rate. 



Figure Bla shows the dust and gas temperature obtained with the thermal balance calculation. As expected, the 
gas temperature is much higher than the dust temperature for low optical depth, while both temperatures agree 
well at high optical depth and density (Doty & Neufeld 1997). The mean of the gas temperature obtained with 
different PDR codes is given by a dotted line, along with one standard deviation given by the gray shaded region. 
The results of our code agrees very well with the results of the PDR codes, while the calculation only takes a few seconds. 

Figure Bib gives the abundances of molecules involved in heating or cooling processes. As Meijerink & Spaans 
(2005), we use a "one-line" approximation for the H2 and CO self-shielding and are thus not able to reproduce the 
exact position of the H/H2 transition. Notice the very low water abundance due to the destruction of water by X-rays 
(Stauber et al. 2006). In Fig. Blc, the rates of each heating process is given along with the total heating rate. At low 
optical depth, photoelectric heating and H2 -pumping are the main heating sources, while dust-gas coupling heats 
the gas at high optical depth (r > 6 in the example). X-ray heating does not contribute signiflcantly to the heating 
rate, but the destruction the important coolant H2O by X-rays slightly enhances the temperature. Fig. Bid shows 
rates of different cooling processes. At the edge of the calculated region, H2 ro-vibrational lines govern the cooling. 
[01] 63 /im and molecular cooling is important at higher optical depth. Due to the low abundance of water, the main 
molecular coolant is CO. 



C. CO+ PRODUCTION IN A MIXING LAYER 

As a toy-model for the mixing layer between the warm and ionized outflow and the envelope, we assume the 
outflow material to be in the form of electrons (e~) and ionized hydrogen (11+ ). The chemical evolution of a parcel 
of gas consisting of a mixture of outflow and inflow material is modeled similarly to Sect. 2. Free parameters are the 
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temperature and the initial abundances of e^ and H+. The peak fractional abundance of the temporal evolution of 
C0+ is given in Fig. CI together with the temporal evolution of the fractional abundance for selected temperatures 
and initial ionizations. 



For temperatures above 60 K, formaldehyde (H2CO) is evaporated from ice mantles. The reaction H2CO + 
H+ C0+ + H2 + H forms C0+ and it is destroyed by the reaction with H2. For very high initial ioniza- 
tion fractions (> 10^^), the dissociative recombination C0+ + e~ ^ C + O with rate coefficient k cx 1/VT 
becomes important and the peak fractional abundance is weakly temperature dependent. Due to the short recombina- 
tion time-scale of e^ with H+ and grains, the C0+ abundance peaks at young chemical ages and then decreases quickly. 



a.) Peak fractional abundance of C0+ b.) Temporal evolution of CO^ 




10 20 50 100 200 500 1000 lo" 10' 10 ' 10° 10" 10"' 

Temperature [K] Time [s] 



Fig. CI. — a.) Peak fractional abundance of the temporal evolution of CO+ for different gas temperatures and initial ionization fractions 
is displayed by marked isocontours and gray scale. The density was chosen to be lO'^ cm""^. b.) Temporal evolution of the fractional 
abundance at points A, B, C, D (black lines). For a temperature of 82 K and an ionization of 10"'' (point B), densities of 10^ cm"'', 10^ 
cm""' and 10^ cm"'' are given in red lines. 



The column density along the mixing layer is A^(CO^) ^ L ■ n ■ Xqq+, with the width of the mixing layer L, the 
gas density n and the peak fractional abundance Xqq+. We set L = va ■ t, with the Alfven velocity va < 10 km s"^ 
and the time t, defined by the point of evolution where the abundance of CO^ has dropped by 1 order of magnitude 
compared to the peak. The peak fractional abundance Xqq+ and the product t ■ n are approximately independent 
of density as Fig. Clb shows. On the other hand, t ■ Xqq+ is roughly independent of temperature and peaks at an 
electron fraction of about 10"^, where < « 5 x 10^ s and a;co+ ~ 2 x 10^^° for a density of 10^ cm~^. We obtain for 
the upper limit of the column denisty due to mixing A^(CO~'') Ki va ■ t ■ n ■ Xqq+ « 10* cm~^. In the FUV model of 
Sect. 3, the column density of C0+ along constant z is of order 10^° cm~^. We conclude that mixing alone cannot 
reproduce the observed amount of C0+. 
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